home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zbesh.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  6.8 KB  |  184 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((hpi 1.5707963267948966))
  12.   (declare (type double-float hpi))
  13.   (defun zbesh (zr zi fnu kode m n cyr cyi nz ierr)
  14.     (declare (type double-float zr zi fnu)
  15.              (type (simple-array double-float (*)) cyr cyi)
  16.              (type f2cl-lib:integer4 kode m n nz ierr))
  17.     (prog ((i 0) (inu 0) (inuh 0) (ir 0) (k 0) (k1 0) (k2 0) (mm 0) (mr 0)
  18.            (nn 0) (nuf 0) (nw 0) (aa 0.0) (alim 0.0) (aln 0.0) (arg 0.0)
  19.            (az 0.0) (dig 0.0) (elim 0.0) (fmm 0.0) (fn 0.0) (fnul 0.0)
  20.            (rhpi 0.0) (rl 0.0) (r1m5 0.0) (sgn 0.0) (str 0.0) (tol 0.0)
  21.            (ufl 0.0) (zni 0.0) (znr 0.0) (zti 0.0) (bb 0.0) (ascle 0.0)
  22.            (rtol 0.0) (atol 0.0) (sti 0.0) (csgnr 0.0) (csgni 0.0))
  23.       (declare
  24.        (type double-float csgni csgnr sti atol rtol ascle bb zti znr zni ufl
  25.         tol str sgn r1m5 rl rhpi fnul fn fmm elim dig az arg aln alim aa)
  26.        (type f2cl-lib:integer4 nw nuf nn mr mm k2 k1 k ir inuh inu i))
  27.       (setf ierr 0)
  28.       (setf nz 0)
  29.       (if (and (= zr 0.0) (= zi 0.0)) (setf ierr 1))
  30.       (if (< fnu 0.0) (setf ierr 1))
  31.       (if (or (< m 1) (> m 2)) (setf ierr 1))
  32.       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  33.       (if (< n 1) (setf ierr 1))
  34.       (if (/= ierr 0) (go end_label))
  35.       (setf nn n)
  36.       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  37.       (setf k1 (f2cl-lib:i1mach 15))
  38.       (setf k2 (f2cl-lib:i1mach 16))
  39.       (setf r1m5 (f2cl-lib:d1mach 5))
  40.       (setf k (f2cl-lib:int (min (abs k1) (abs k2))))
  41.       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  42.       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  43.       (setf aa (* r1m5 k1))
  44.       (setf dig (min aa 18.0))
  45.       (setf aa (* aa 2.303))
  46.       (setf alim (+ elim (max (- aa) -41.45)))
  47.       (setf fnul (+ 10.0 (* 6.0 (- dig 3.0))))
  48.       (setf rl (+ (* 1.2 dig) 3.0))
  49.       (setf fn (+ fnu (f2cl-lib:int-sub nn 1)))
  50.       (setf mm (f2cl-lib:int-sub 3 m m))
  51.       (setf fmm (coerce (the f2cl-lib:integer4 mm) 'double-float))
  52.       (setf znr (* fmm zi))
  53.       (setf zni (* (- fmm) zr))
  54.       (setf az (zabs zr zi))
  55.       (setf aa (/ 0.5 tol))
  56.       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
  57.       (setf aa (min aa bb))
  58.       (f2cl-lib:fformat t (("~A~%")) "az, aa = " az aa)
  59.       (f2cl-lib:fformat t (("~A~%")) "fn     = " fn)
  60.       (if (> az aa) (go label260))
  61.       (if (> fn aa) (go label260))
  62.       (setf aa (f2cl-lib:fsqrt aa))
  63.       (if (> az aa) (setf ierr 3))
  64.       (if (> fn aa) (setf ierr 3))
  65.       (setf ufl (* (f2cl-lib:d1mach 1) 1000.0))
  66.       (if (< az ufl) (go label230))
  67.       (if (> fnu fnul) (go label90))
  68.       (if (<= fn 1.0) (go label70))
  69.       (if (> fn 2.0) (go label60))
  70.       (if (> az tol) (go label70))
  71.       (setf arg (* 0.5 az))
  72.       (setf aln (* (- fn) (f2cl-lib:flog arg)))
  73.       (if (> aln elim) (go label230))
  74.       (go label70)
  75.      label60
  76.       (multiple-value-bind
  77.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  78.            var-11)
  79.           (zuoik znr zni fnu kode 2 nn cyr cyi nuf tol elim alim)
  80.         (declare
  81.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9 var-10
  82.           var-11))
  83.         (setf nuf var-8))
  84.       (if (< nuf 0) (go label230))
  85.       (setf nz (f2cl-lib:int-add nz nuf))
  86.       (setf nn (f2cl-lib:int-sub nn nuf))
  87.       (if (= nn 0) (go label140))
  88.      label70
  89.       (if (or (< znr 0.0) (and (= znr 0.0) (< zni 0.0) (= m 2))) (go label80))
  90.       (multiple-value-bind
  91.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10)
  92.           (zbknu znr zni fnu kode nn cyr cyi nz tol elim alim)
  93.         (declare
  94.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9 var-10))
  95.         (setf nz var-7))
  96.       (go label110)
  97.      label80
  98.       (setf mr (f2cl-lib:int-sub mm))
  99.       (multiple-value-bind
  100.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  101.            var-11 var-12 var-13)
  102.           (zacon znr zni fnu kode mr nn cyr cyi nw rl fnul tol elim alim)
  103.         (declare
  104.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9 var-10
  105.           var-11 var-12 var-13))
  106.         (setf nw var-8))
  107.       (if (< nw 0) (go label240))
  108.       (setf nz nw)
  109.       (go label110)
  110.      label90
  111.       (setf mr 0)
  112.       (if (and (>= znr 0.0) (or (/= znr 0.0) (>= zni 0.0) (/= m 2)))
  113.           (go label100))
  114.       (setf mr (f2cl-lib:int-sub mm))
  115.       (if (or (/= znr 0.0) (>= zni 0.0)) (go label100))
  116.       (setf znr (- znr))
  117.       (setf zni (- zni))
  118.      label100
  119.       (multiple-value-bind
  120.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  121.            var-11)
  122.           (zbunk znr zni fnu kode mr nn cyr cyi nw tol elim alim)
  123.         (declare
  124.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9 var-10
  125.           var-11))
  126.         (setf nw var-8))
  127.       (if (< nw 0) (go label240))
  128.       (setf nz (f2cl-lib:int-add nz nw))
  129.      label110
  130.       (setf sgn (coerce (f2cl-lib:dsign hpi (- fmm)) 'double-float))
  131.       (setf inu (f2cl-lib:int fnu))
  132.       (setf inuh (the f2cl-lib:integer4 (truncate inu 2)))
  133.       (setf ir (f2cl-lib:int-sub inu (f2cl-lib:int-mul 2 inuh)))
  134.       (setf arg (* (- fnu (f2cl-lib:int-sub inu ir)) sgn))
  135.       (setf rhpi (/ 1.0 sgn))
  136.       (setf csgni (* rhpi (cos arg)))
  137.       (setf csgnr (* (- rhpi) (sin arg)))
  138.       (if (= (mod inuh 2) 0) (go label120))
  139.       (setf csgnr (- csgnr))
  140.       (setf csgni (- csgni))
  141.      label120
  142.       (setf zti (- fmm))
  143.       (setf rtol (/ 1.0 tol))
  144.       (setf ascle (* ufl rtol))
  145.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  146.                     ((> i nn) nil)
  147.         (tagbody
  148.           (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
  149.           (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
  150.           (setf atol 1.0)
  151.           (if (> (max (abs aa) (abs bb)) ascle) (go label135))
  152.           (setf aa (* aa rtol))
  153.           (setf bb (* bb rtol))
  154.           (setf atol tol)
  155.          label135
  156.           (setf str (- (* aa csgnr) (* bb csgni)))
  157.           (setf sti (+ (* aa csgni) (* bb csgnr)))
  158.           (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 n))) (* str atol))
  159.           (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 n))) (* sti atol))
  160.           (setf str (* (- csgni) zti))
  161.           (setf csgni (* csgnr zti))
  162.           (setf csgnr str)
  163.          label130))
  164.       (go end_label)
  165.      label140
  166.       (if (< znr 0.0) (go label230))
  167.       (go end_label)
  168.      label230
  169.       (setf nz 0)
  170.       (setf ierr 2)
  171.       (go end_label)
  172.      label240
  173.       (if (= nw -1) (go label230))
  174.       (setf nz 0)
  175.       (setf ierr 5)
  176.       (go end_label)
  177.      label260
  178.       (setf nz 0)
  179.       (setf ierr 4)
  180.       (go end_label)
  181.      end_label
  182.       (return (values nil nil nil nil nil nil nil nil nz ierr)))))
  183.  
  184.